/*==============================================================================
    FAC MONETARY POLICY ROBUSTNESS TESTS
    
    Three state-of-the-art robustness methods:
    1. Local Projections (Jordà 2005)
    2. Lewbel (2012) Heteroskedasticity-Based IV
    3. Difference-in-Differences (Regime Breaks)
    
    Data: fac_master_FINAL.csv (231 quarters, 1968 Q1 - 2025 Q3)
    
    Author: Ciarán
    Date: December 2024
==============================================================================*/

clear all
set more off
set scheme s1mono

* Set working directory - adjust path as needed
cd "/Users/ciaran/PycharmProjects/fac"

* Install required packages
cap ssc install outreg2
cap ssc install ivreg2h
cap ssc install ivreg2
cap ssc install ranktest
cap ssc install freduse

* Create output directories
cap mkdir "results"
cap mkdir "results/local_projections"
cap mkdir "results/lewbel_iv"
cap mkdir "results/did"
cap mkdir "results/figures"


/*------------------------------------------------------------------------------
    SECTION 1: DATA PREPARATION
------------------------------------------------------------------------------*/

* Import FAC master data
import delimited "fac_master_FINAL.csv", clear varnames(1)

* Create quarterly time variable
gen yq = yq(year, quarter)
format yq %tq
tsset yq

* Rename score to fac_score for clarity
rename score fac_score

* Create lagged FAC scores
gen L1_fac = L.fac_score
gen L2_fac = L2.fac_score
gen L3_fac = L3.fac_score
gen L4_fac = L4.fac_score

* Save FAC data
tempfile facdata
save `facdata'

/*------------------------------------------------------------------------------
    FRED DATA DOWNLOAD
    
    Alternative 1: Use freduse (requires internet)
    Alternative 2: Use existing fac_fred_merged.dta
------------------------------------------------------------------------------*/

* Option A: Download from FRED (uncomment if needed)
/*
freduse FEDFUNDS CPIAUCSL UNRATE GDPC1, clear

* Process Fed Funds
rename FEDFUNDS fedfunds
gen yq = qofd(daten)
format yq %tq
collapse (mean) fedfunds, by(yq)
tempfile ff
save `ff'

* Process CPI → Inflation
freduse CPIAUCSL, clear
gen yq = qofd(daten)
format yq %tq
collapse (mean) CPIAUCSL, by(yq)
tsset yq
gen inflation = 400 * (CPIAUCSL / L.CPIAUCSL - 1)
tempfile cpi
save `cpi'

* Process Unemployment
freduse UNRATE, clear
gen yq = qofd(daten)
format yq %tq
collapse (mean) unemployment=UNRATE, by(yq)
tempfile ur
save `ur'

* Merge FRED data
use `ff', clear
merge 1:1 yq using `cpi', nogen
merge 1:1 yq using `ur', nogen
tempfile freddata
save `freddata'
*/

* Option B: Use pre-existing FRED data (assumes fred_macro.dta exists)
* If you have fac_fred_merged.dta from previous work:
cap confirm file "fred_macro.dta"
if _rc != 0 {
    di as error "FRED data file not found. Creating from FRED API..."
    
    * Download Federal Funds Rate
    freduse FEDFUNDS, clear
    gen yq = qofd(daten)
    format yq %tq
    collapse (mean) fedfunds=FEDFUNDS, by(yq)
    tempfile ff
    save `ff'
    
    * Download CPI for inflation
    freduse CPIAUCSL, clear
    gen yq = qofd(daten)
    format yq %tq
    collapse (mean) cpi=CPIAUCSL, by(yq)
    tsset yq
    gen inflation = 400 * (cpi / L.cpi - 1)
    drop cpi
    tempfile inf
    save `inf'
    
    * Download Unemployment
    freduse UNRATE, clear
    gen yq = qofd(daten)
    format yq %tq
    collapse (mean) unemployment=UNRATE, by(yq)
    tempfile ur
    save `ur'
    
    * Download GDP for output gap
    freduse GDPC1, clear
    gen yq = qofd(daten)
    format yq %tq
    collapse (mean) gdp=GDPC1, by(yq)
    tsset yq
    gen gdp_growth = 400 * (gdp / L.gdp - 1)
    drop gdp
    tempfile gdp
    save `gdp'
    
    * Merge all FRED series
    use `ff', clear
    merge 1:1 yq using `inf', nogen
    merge 1:1 yq using `ur', nogen
    merge 1:1 yq using `gdp', nogen
    
    save "fred_macro.dta", replace
}
else {
    use "fred_macro.dta", clear
}

* Merge with FAC data
merge 1:1 yq using `facdata', keep(match) nogen

* Create dependent variables
tsset yq
gen dfedfunds = D.fedfunds
gen F1_fedfunds = F.fedfunds
gen F2_fedfunds = F2.fedfunds
gen F4_fedfunds = F4.fedfunds
gen F8_fedfunds = F8.fedfunds
gen F12_fedfunds = F12.fedfunds

* Create change in Fed Funds at various horizons for LP
forvalues h = 0/12 {
    gen dff_h`h' = F`h'.fedfunds - L.fedfunds
}

* Create regime dummy variables
gen post_volcker = (year >= 1980)
gen post_glb = (year >= 2000)
gen post_gfc = (year >= 2010)

* Create interaction terms
gen fac_post_volcker = L1_fac * post_volcker
gen fac_post_glb = L1_fac * post_glb
gen fac_post_gfc = L1_fac * post_gfc

/*------------------------------------------------------------------------------
    VARIABLE LABELS
------------------------------------------------------------------------------*/

label variable fac_score "FAC Monetary Policy Score"
label variable L1_fac "FAC Score (t-1)"
label variable L2_fac "FAC Score (t-2)"
label variable L3_fac "FAC Score (t-3)"
label variable L4_fac "FAC Score (t-4)"
label variable fedfunds "Federal Funds Rate (%)"
label variable dfedfunds "Δ Federal Funds Rate"
label variable inflation "Inflation Rate (%)"
label variable unemployment "Unemployment Rate (%)"
label variable gdp_growth "GDP Growth (%)"
label variable post_volcker "Post-Volcker (≥1980)"
label variable post_glb "Post-GLB (≥2000)"
label variable post_gfc "Post-GFC (≥2010)"
label variable fac_post_volcker "FAC Score × Post-Volcker"
label variable fac_post_glb "FAC Score × Post-GLB"
label variable fac_post_gfc "FAC Score × Post-GFC"

forvalues h = 0/12 {
    label variable dff_h`h' "Δ Fed Funds (h=`h')"
}

* Save analysis dataset
save "fac_analysis.dta", replace


/*==============================================================================
    ROBUSTNESS TEST 1: LOCAL PROJECTIONS (Jordà 2005)
    
    IRF without VAR assumptions
    Newey-West standard errors for serial correlation
==============================================================================*/

di _n(3) "======================================================================"
di "LOCAL PROJECTIONS: FAC Score → Federal Funds Rate"
di "======================================================================"

* Storage matrices for IRF
matrix IRF = J(13, 5, .)
matrix colnames IRF = v1 v2 v3 v4 v5

* Run local projections at each horizon
forvalues h = 0/12 {
    
    * Newey-West with lag = horizon + 1
    local nw_lag = `h' + 1
    
    * LP regression
    newey dff_h`h' L1_fac L.inflation L.unemployment, lag(`nw_lag')
    
    * Store results
    matrix IRF[`h'+1, 1] = `h'
    matrix IRF[`h'+1, 2] = _b[L1_fac]
    matrix IRF[`h'+1, 3] = _se[L1_fac]
    matrix IRF[`h'+1, 4] = _b[L1_fac] - 1.96 * _se[L1_fac]
    matrix IRF[`h'+1, 5] = _b[L1_fac] + 1.96 * _se[L1_fac]
    
    * Store for outreg2 at key horizons
    if inlist(`h', 0, 1, 4, 8, 12) {
        est store lp_h`h'
    }
}

* Create IRF dataset for plotting
clear
svmat IRF

* Rename for clarity
rename IRF1 horizon
rename IRF2 beta
rename IRF3 se_nw
rename IRF4 ci_lo
rename IRF5 ci_hi

* Generate plot
gen zero = 0

twoway (rarea ci_lo ci_hi horizon, color(gs14) lwidth(none)) ///
       (line beta horizon, lcolor(navy) lwidth(medium)) ///
       (line zero horizon, lpattern(dash) lcolor(gs8)), ///
       ytitle("Response of Fed Funds Rate (pp)", size(medium)) ///
       xtitle("Quarters After FAC Meeting", size(medium)) ///
       title("Impulse Response: FAC Score → Fed Funds", size(large)) ///
       subtitle("Local Projections with Newey-West SEs", size(small)) ///
       legend(off) ///
       xlabel(0(2)12) ///
       ylabel(, angle(horizontal)) ///
       note("Shaded area: 95% confidence interval" ///
            "Controls: lagged inflation and unemployment")
            
graph export "results/figures/irf_local_projections.png", replace width(1200)

* Export IRF data
export delimited using "results/local_projections/irf_data.csv", replace

* Reload main data for regression output
use "fac_analysis.dta", clear

* Re-run and export table for key horizons
forvalues h = 0/12 {
    local nw_lag = `h' + 1
    newey dff_h`h' L1_fac L.inflation L.unemployment, lag(`nw_lag')
    
    if `h' == 0 {
        outreg2 using "results/local_projections/table_local_projections.doc", ///
            replace word label dec(3) ///
            title("Table: Local Projections - FAC Score → Federal Funds Rate") ///
            ctitle("h=`h'") ///
            addtext(Controls, "Inflation, Unemployment", SE Type, "Newey-West")
    }
    else if inlist(`h', 1, 4, 8, 12) {
        outreg2 using "results/local_projections/table_local_projections.doc", ///
            append word label dec(3) ///
            ctitle("h=`h'")
    }
}

di _n "Local Projections complete. Results saved to results/local_projections/"


/*==============================================================================
    ROBUSTNESS TEST 2: LEWBEL (2012) HETEROSKEDASTICITY-BASED IV
    
    Internal instruments from heteroskedasticity
    Tests endogeneity without external instruments
==============================================================================*/

di _n(3) "======================================================================"
di "LEWBEL IV: Heteroskedasticity-Based Identification"
di "======================================================================"

* Standard OLS for comparison
reg dfedfunds L1_fac L.inflation L.unemployment, robust
est store ols_base

outreg2 using "results/lewbel_iv/table_lewbel_comparison.doc", ///
    replace word label dec(3) ///
    title("Table: OLS vs Lewbel IV Estimates") ///
    ctitle("OLS") ///
    addtext(Method, "OLS", Instruments, "None", Period, "1968-2025")

* Lewbel IV estimation
* ivreg2h treats L1_fac as potentially endogenous
* Uses heteroskedasticity in first-stage residuals as instruments

cap ivreg2h dfedfunds L.inflation L.unemployment (L1_fac = ), robust first

if _rc == 0 {
    * Store Lewbel results
    est store lewbel_iv
    
    * Manual extraction since ivreg2h doesn't directly support outreg2
    local b_fac = _b[L1_fac]
    local se_fac = _se[L1_fac]
    local n = e(N)
    
    * Re-run and append (workaround for outreg2 compatibility)
    * Use ivreg2 with generated instruments as approximation
    
    * Generate heteroskedasticity-based instruments manually
    qui reg L1_fac L.inflation L.unemployment
    predict resid_first, resid
    
    * Create Lewbel-style instruments: (Z - mean(Z)) * resid
    foreach var in L1inflation L1unemployment {
        cap gen `var' = L.inflation if "`var'" == "L1inflation"
        cap replace `var' = L.unemployment if "`var'" == "L1unemployment"
    }
    
    egen mean_inf = mean(L.inflation)
    egen mean_unemp = mean(L.unemployment)
    
    gen z1 = (L.inflation - mean_inf) * resid_first
    gen z2 = (L.unemployment - mean_unemp) * resid_first
    
    * IV regression with generated instruments
    ivreg2 dfedfunds L.inflation L.unemployment (L1_fac = z1 z2), robust first
    est store iv_manual
    
    outreg2 using "results/lewbel_iv/table_lewbel_comparison.doc", ///
        append word label dec(3) ///
        ctitle("Lewbel IV") ///
        addtext(Method, "Lewbel IV", Instruments, "Heteroskedasticity", Period, "1968-2025")
    
    * First stage results
    di _n "First Stage Results:"
    reg L1_fac z1 z2 L.inflation L.unemployment, robust
    test z1 z2
    local fstat = r(F)
    di "First-stage F-statistic: " %5.2f `fstat'
    
    * Save first stage
    outreg2 using "results/lewbel_iv/table_first_stage.doc", ///
        replace word label dec(3) ///
        title("Table: Lewbel IV First Stage") ///
        ctitle("FAC Score (t-1)") ///
        addstat(F-statistic, `fstat') ///
        addtext(Period, "1968-2025")
}
else {
    di as error "ivreg2h failed. Running manual Lewbel procedure..."
    
    * Manual Lewbel implementation
    qui reg L1_fac L.inflation L.unemployment
    predict resid_first, resid
    
    egen mean_inf = mean(L.inflation)
    egen mean_unemp = mean(L.unemployment)
    
    gen z1 = (L.inflation - mean_inf) * resid_first
    gen z2 = (L.unemployment - mean_unemp) * resid_first
    
    * 2SLS with generated instruments
    ivregress 2sls dfedfunds L.inflation L.unemployment (L1_fac = z1 z2), robust first
    
    outreg2 using "results/lewbel_iv/table_lewbel_comparison.doc", ///
        append word label dec(3) ///
        ctitle("Lewbel IV") ///
        addtext(Method, "Manual Lewbel", Instruments, "Heteroskedasticity", Period, "1968-2025")
}

* Endogeneity test (Durbin-Wu-Hausman via ivreg2)
di _n "Endogeneity Test:"
qui ivreg2 dfedfunds L.inflation L.unemployment (L1_fac = z1 z2), robust endog(L1_fac)
di "Endogeneity test p-value: " %5.4f e(estatp)

di _n "Lewbel IV complete. Results saved to results/lewbel_iv/"


/*==============================================================================
    ROBUSTNESS TEST 3: DIFFERENCE-IN-DIFFERENCES (Regime Breaks)
    
    Three structural breaks:
    1. Volcker (1980) - Monetary policy revolution
    2. GLB (2000) - Gramm-Leach-Bliley deregulation
    3. GFC (2010) - Post-crisis regime
==============================================================================*/

di _n(3) "======================================================================"
di "DIFFERENCE-IN-DIFFERENCES: Regime Breaks"
di "======================================================================"

*------------------------------------------------------------------------------
* 3A. Volcker Break (1980)
*------------------------------------------------------------------------------

di _n "--- Volcker Break (1980) ---"

* Baseline: FAC effect on ΔFedfunds
reg dfedfunds L1_fac L.inflation L.unemployment, robust
est store did_base

outreg2 using "results/did/table_did_regimes.doc", ///
    replace word label dec(3) ///
    title("Table: Difference-in-Differences - Regime Breaks") ///
    ctitle("Baseline") ///
    addtext(Break, "None", Controls, "Yes", Period, "1968-2025")

* Volcker interaction
reg dfedfunds c.L1_fac##c.post_volcker L.inflation L.unemployment, robust
est store did_volcker

* Marginal effects
margins, dydx(L1_fac) at(post_volcker=(0 1))
marginsplot, ///
    title("FAC Effect: Pre vs Post-Volcker") ///
    ytitle("Effect on Δ Fed Funds Rate") ///
    xtitle("") ///
    xlabel(0 "Pre-1980" 1 "Post-1980") ///
    recast(bar) ///
    plotopts(barwidth(0.5) color(navy%70)) ///
    ciopts(lcolor(black))
graph export "results/figures/margins_volcker.png", replace width(1200)

outreg2 using "results/did/table_did_regimes.doc", ///
    append word label dec(3) ///
    ctitle("Volcker") ///
    addtext(Break, "1980", Controls, "Yes", Period, "1968-2025")

*------------------------------------------------------------------------------
* 3B. GLB Break (2000)
*------------------------------------------------------------------------------

di _n "--- GLB Break (2000) ---"

reg dfedfunds c.L1_fac##c.post_glb L.inflation L.unemployment, robust
est store did_glb

margins, dydx(L1_fac) at(post_glb=(0 1))
marginsplot, ///
    title("FAC Effect: Pre vs Post-GLB") ///
    ytitle("Effect on Δ Fed Funds Rate") ///
    xtitle("") ///
    xlabel(0 "Pre-2000" 1 "Post-2000") ///
    recast(bar) ///
    plotopts(barwidth(0.5) color(navy%70)) ///
    ciopts(lcolor(black))
graph export "results/figures/margins_glb.png", replace width(1200)

outreg2 using "results/did/table_did_regimes.doc", ///
    append word label dec(3) ///
    ctitle("GLB") ///
    addtext(Break, "2000", Controls, "Yes", Period, "1968-2025")

*------------------------------------------------------------------------------
* 3C. GFC Break (2010)
*------------------------------------------------------------------------------

di _n "--- GFC Break (2010) ---"

reg dfedfunds c.L1_fac##c.post_gfc L.inflation L.unemployment, robust
est store did_gfc

margins, dydx(L1_fac) at(post_gfc=(0 1))
marginsplot, ///
    title("FAC Effect: Pre vs Post-GFC") ///
    ytitle("Effect on Δ Fed Funds Rate") ///
    xtitle("") ///
    xlabel(0 "Pre-2010" 1 "Post-2010") ///
    recast(bar) ///
    plotopts(barwidth(0.5) color(navy%70)) ///
    ciopts(lcolor(black))
graph export "results/figures/margins_gfc.png", replace width(1200)

outreg2 using "results/did/table_did_regimes.doc", ///
    append word label dec(3) ///
    ctitle("GFC") ///
    addtext(Break, "2010", Controls, "Yes", Period, "1968-2025")

*------------------------------------------------------------------------------
* 3D. All Breaks Combined
*------------------------------------------------------------------------------

di _n "--- All Breaks Combined ---"

reg dfedfunds c.L1_fac##c.post_volcker c.L1_fac##c.post_glb c.L1_fac##c.post_gfc ///
    L.inflation L.unemployment, robust
est store did_all

outreg2 using "results/did/table_did_regimes.doc", ///
    append word label dec(3) ///
    ctitle("All Breaks") ///
    addtext(Breaks, "1980, 2000, 2010", Controls, "Yes", Period, "1968-2025")

*------------------------------------------------------------------------------
* Summary Margins Plot (all periods)
*------------------------------------------------------------------------------

* Create period variable for comprehensive margins
gen period = 1 if year < 1980
replace period = 2 if year >= 1980 & year < 2000
replace period = 3 if year >= 2000 & year < 2010
replace period = 4 if year >= 2010

label define period_lbl 1 "1968-1979" 2 "1980-1999" 3 "2000-2009" 4 "2010-2025"
label values period period_lbl

reg dfedfunds c.L1_fac##i.period L.inflation L.unemployment, robust
margins, dydx(L1_fac) at(period=(1 2 3 4))
marginsplot, ///
    title("FAC Effect Across Monetary Policy Regimes") ///
    ytitle("Effect on Δ Fed Funds Rate") ///
    xtitle("") ///
    xlabel(1 "Pre-Volcker" 2 "Volcker-GLB" 3 "GLB-GFC" 4 "Post-GFC", angle(15)) ///
    recast(bar) ///
    plotopts(barwidth(0.5) color(navy%70)) ///
    ciopts(lcolor(black)) ///
    note("Bars show marginal effect of FAC Score (t-1) on Δ Fed Funds" ///
         "Error bars: 95% confidence intervals")
graph export "results/figures/margins_all_periods.png", replace width(1200)

di _n "Difference-in-Differences complete. Results saved to results/did/"


/*==============================================================================
    SUMMARY OUTPUT
==============================================================================*/

di _n(3) "======================================================================"
di "ROBUSTNESS TESTS COMPLETE"
di "======================================================================"

di _n "Output files created:"
di "  results/local_projections/table_local_projections.doc"
di "  results/local_projections/irf_data.csv"
di "  results/figures/irf_local_projections.png"
di ""
di "  results/lewbel_iv/table_lewbel_comparison.doc"
di "  results/lewbel_iv/table_first_stage.doc"
di ""
di "  results/did/table_did_regimes.doc"
di "  results/figures/margins_volcker.png"
di "  results/figures/margins_glb.png"
di "  results/figures/margins_gfc.png"
di "  results/figures/margins_all_periods.png"

* Save final dataset
save "fac_analysis_final.dta", replace

di _n "Analysis dataset saved: fac_analysis_final.dta"
di "======================================================================"

log close, replace

/*==============================================================================
    END OF DO-FILE
==============================================================================*/
